Fluctuations, line tensions, and correlation times of nanoscale islands on surfaces 
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We analyze in detail the fluctuations and correlations of the (spatial) Fourier modes of nano-scale 
single-layer islands on (111) fee crystal surfaces. We analytically show that the Fourier modes of 
the fluctuations couple due to the anisotropy of the crystal, changing the power spectrum of the 
fluctuations, and that the actual eigenmodes of the fluctuations are the appropriate linear combina- 
tions of the Fourier modes. Using kinetic Monte Carlo simulations with bond-counting parameters 
that best match realistic energy barriers for hopping rates, we deduce absolute line tensions as a 
function of azimuthal orientation from the analyses of the fluctuation of each individual mode. The 
autocorrelation functions of these modes give the scaling of the correlation times with wavelength, 
providing us with the rate-limiting kinetics driving the fluctuations, here step-edge diffusion. The 
results for the energetic parameters are in reasonable agreement with available experimental data 
for Pb(lll) surfaces, and we compare the correlation times of island-edge fluctuations to relaxation 
times of quenched Pb crystallites. 
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I. INTRODUCTION 

Nanoscale islands consisting of 10 2 - 10 5 atoms have 
captured great interest over recent years for a variety 
of reasons. From a practical standpoint, they provide 
a precursor to the formation of quantum dots, which, 
if assembled in a controlled way, can serve as the basic 
ingredients of nano-scale electronic and mechanical de- 
vices. Many crystallites or nanomounds are best viewed 
as "wedding-cake" -like stacks of such islandsi They are 
the intermediary between a flat surface and a small three- 
dimensional structure. In contrast to steps, which re- 
quire vicinal surfaces^ that often must be well character- 
ized over mcsoscopic regions, islands can be studied in 
smaller-scale regions that are flat only locally. 

Of particular interest to us are the shape and the fluc- 
tuations of the perimeter of these islands. The shape 
provides information about the line tension or step free 
energy per length, from which one can compute the step 
stiffness that describes the "inertial" properties of steps. 
The "dipole" mode of these fluctuations are known to 
underlie the diffusion of such islands, a concept now used 
routinely in simulations^^ However, shorter-wavelength 
modes are also of great interest, since they can be cor- 
related with similar fluctuations of steps and provide a 
way to assess, again, the stiffness of the step and also the 
kinetic or atomistic diffusion coefficient associated with 
the mechanism that dominates the atomistic processes 
underlying the fluctuations. Until recently, attention was 
limited to structures for which crystal anisotropy could 
seemingly be ignored. 

Here we pay particular attention to the role of the 
inevitable anisotropy of crystal surfaces, which around 
room temperature or even above it is typically sufficiently 
strong that it should apparently be taken into account in 
order to correctly characterize the morphology of the var- 
ious (near) equilibrium structures appearing on surfaces 



and their dynamics. In this paper we focus on the line 
tension and stiffness and their orientation dependence; 
we give an analytic method to calculating these physical 
parameters from the fluctuation of nanoscale islands. 

The little experimental data on such systems involve 
runs of worrisome duration or use probes that provide 
scanned rather than instantaneous images. To gener- 
ate fully-characterized data, we turned to kinetic Monte 
Carlo (KMC) simulations to mimic the equilibrium fluc- 
tuations of islands. These simulations are the input of 
our analytic theory which, starting from the excess free 
energy corresponding to the capillary wave fluctuations 
of the island edge, provides the eigenmodes of these fluc- 
tuations. Since the 2D Wulff plot relating the equilibrium 
island shape and the line tension in the azimuthal direc- 
tions on the surface provides only relative line tensions for 
various orientations, a key problem is always the deter- 
mination of the chemical potential A of the island edge, 
which then produces an absolute relation. This poten- 
tial can be determined with surprisingly good (~ 10%) 
accuracy from the spectrum of the modes of the system. 
We compare these eigenmodes and the simple Fourier 
modes of the fluctuations and reach the (perhaps) sur- 
prising conclusion that the anisotropy only affects the 
longer wavelength modes. 

Another aim of the paper is to examine the correlation 
of the fluctuations of the Fourier modes and thereby to 
find the rate-limiting process driving the fluctuations in 
a fairly realistic model. For our KMC simulations we 
sought a system for which one could compute hop rates 
with good accuracy and for which there was quantitative 
experimental data with which to compare. Accordingly, 
we have chosen Pb(lll) so as to be able to compare with 
intriguing recent experiments by Thurmer et al& This 
analysis gives the scaling of the correlation time with the 
wavelength, that is the dynamic exponent z, and provides 
us with characteristic times measured not only in MC 
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FIG. 1: (Color online) Geometry of the MC simulation. The 
approximate mean radius R of the island and the radius R c 
of the container are illustrated. 

steps, but in real time. Thus, we can compare directly 
with experiments and extrapolate to different structures 
from the simple one considered here. 

Utilizing direct surface imaging techniques, especially 
scanning tunneling microscopy (STM), several attempts 
have been made to measure and calculate step energies. 
From a theoretical viewpoint the various methods that 
used the experimental data for calculations can be broken 
down into two main groups. The first is based on a lattice 
model which relates the island shape (radius and curva- 
ture) to the temperature dependence of the free energy 
and stiffness of the Ising model in the low-temperature 
expansion, usually in high symmetry directions. By fit- 
ting the functional shape of the free energy with vary- 
ing temperature on the experimental data determined by 
the equilibrium island shaped gives the Ising kink energy, 
which in turn provides the step energies and stiffnesses. 
However, limitations of the Ising model to describe sur- 
face structure have recently been notedi 

The other method is based on a step continuum 
model which makes use of stochastic differential equa- 
tions to describe the fluctuations of straight steps^ or 
island edges&ifi viewed as nearly circular closed-loop 
steps. Thus, the initial calculations for island fluctu- 
ations assumed isotropytii the power spectrum of the 
Fourier modes of the step fluctuations were calculated 
and adapted with appropriate modifications to nearly cir- 
cular island shapes^ If the anisotropy turns out to be 
strong, it cannot be handled as a perturbation; a com- 
plete anisotropic calculation without any such assump- 
tions becomes necessary. 

This challenge was recently taken up by Khare et al.^ 
who give an approximate form for the free energy func- 
tional and calculate the chemical potential integrating 
all the Fourier modes in the system by using the gener- 
alized cquipartition theorem where the modes are buried 



in a sum. However, these modes are coupled, so any one 
mode missing (e.g. due to lack of experimental resolution) 
in the sum can contribute to a deviation from the precise 
value of the chemical potential by itself and through its 
coupling to the other modes as well. In contrast, our ap- 
proach of analyzing individual modes gives more insight 
into the extent to which this coupling should be taken 
into account, and provides the chemical potential in a 
(mathematically) controlled way. 

The autocorrelation function of fluctuations of step 
edges and correlation times have been analyzed theoreti- 
cally in Fourier space based on Langevin formalism^i*« 
and experimentally in the context of straight-step fluctu- 
ations on Si(lll)i*i and Si(001)£ surfaces for relatively 
long wavelengths. The rate-limiting kinetics driving 
these fluctuations arc determined by the dynamic expo- 
nent, which also sets the universality class to which the 
system belongs^ (The roughness exponent is believed 
to be 1/2 in our cases.) The correlation times are the- 
oretically identical to the relaxation times (or, in some 
cases, decay times) of surface features; 17 i 18 such as de- 
cay and near-equilibrium build-up of bulges (of either 
sign) along the step edge that also have (wave)length L. 
Three-dimensional features like mcsoscopic (or smaller) 
wires on surfaces as well as the surface corrugations in 
earlier studies by Mullinsji2i2& are typically described by 
1+1 dimensional models but may involve different, more 
complicated mechanisms driving their fluctuation or de- 
cay. We will compare and discuss these various relaxation 
times in the paper. 

The paper is organized as follows: In the next sec- 
tion we give an analytic solution to the decoupling of the 
Fourier modes of the system into the actual eigenmodes 
and recalculate the free energy functional of the edge fluc- 
tuations. The results and conclusions can be understood 
without the reader's going through this algebra; only the 
result expressed in Eq. fTDjl is used later. In Sec. Ill we 
introduce the KMC simulation and in Sec. IV use its re- 
sults to calculate the chemical potential and line tension. 
In Sec. V we calculate the correlation functions of the 
Fourier modes and deduce the scaling of the correlation 
time with length, the dynamic exponent z. We compare 
with available experimental data for Pb(lll). Sec. VI 
concludes the paper. 



II. FOURIER MODES, EIGENMODES 

The relationship between the equilibrium crystal shape 
and the surface tension or, in our 2D case, between the 
equilibrium island shape and the line tension of its edge 
can be established by the minimization of the free energy 
functional of the island edge. The orientation-dependent 
line tension /3(n) is defined as the work per unit length 
necessary to create the ds line element with normal n 
to the perimeter. The free energy is the integral of this 
work along the whole perimeter. The equilibrium island 
shape at a constant temperature T, number of particles 
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FIG. 2: (Color online) Schematic showing variables used to 
analyze equilibrium island shape 



N, and area £, is determined by the minimization of the 
free energy functional with respect to the shape with the 
constraint that the island area is constant, typically using 
the method of Lagrange multipliers i21 



F[R,R,0] = I (3(n)ds-X j da 
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/3(V(6»)) (R + R ) dB-Xj —de. (1) 



Here the second line is in polar coordinates with e 
the polar angle and R(e) the radius of the equilibrium 
shape. The dot denotes the differentiation with respect 
to the angle, A is the Lagrange multiplier (which actu- 
ally turns out to be the chemical potential), and ip is 
the angle which characterizes the vector normal to the 
shape (see Fig- EJ - ds and da are the line element and 
surface element, respectively. Formally minimizing the 
F = F[R, R, 6] functional, the Euler-Lagrange equation 
gives a relation between the equilibrium island shape 
R(e) and the orientation dependent line tension, 
and between the two angles involved: %j) which depends 
on the polar angle and the equilibrium shape, 1 - 2 - 2 



6F_ 

Jr 



= A 



(r 2 +r 2 Y / - 
arctan 



(2) 



However, in this procedure A is a prefactor and cannot 
be determined, leaving the relation relative. Eq. J2J) is 
the seminal Wulff construction in polar coordinates. 

In order to determine the chemical potential, the ther- 
mal fluctuations of the island edge can be utilized. In this 
case the free energy of the island changes as its shape 
changes due to the fluctuations, and the free energy is 
certainly not at its minimum but depends on the island's 
instantaneous shape. Then the free energy of this in- 
stantaneous shape is the integral over the line elements 



of the shape with their corresponding line tension, which 
changes with time as the orientation of the shape element 
changes: 

L 

'fi(i/){9)) (OR + r? + (-R + rf) '^de (3) 



Here r and r are time dependent and describe the devia- 
tion of the instantaneous shape from the equilibrium one 
as shown in Fig. [3] The angle ip is also time-dependent 
since now it depends not only on R and R as in Eq. J5J), 
but also on r and r. Considering only small deformations 
from the equilibrium shape (as it is usually assumed in 
the capillary wave theory) and also small slope deviations 
from the equilibrium slope R, so that r, r <C R, the Tay- 
lor expansion (both in (3 and in the square root) in these 
small parameters leads to the functional 



F[r,r,0;t]=X 



1 



Rr - Rr 



2R 2 + 2R 2 - RR 



-de. 



(4) 



This functional contains 3 quadratic terms A(9)r 2 , 
Q(9)rr and B(9)r 2 . The cross term Q drops out after 
taking the ensemble average; both the two other terms 
are determined by properties of the equilibrium island 
shape: 



A{9) 
B{9) = 



R 2 



2 R 2 + 2R 2 - RR 

1 R 2 

2 R? + 2R 2 - RR ' 



(5) 
(G) 



and provide the weightings of the fluctuations of the 
deformations characterized by r 2 and f 2 , respectively. 
These deformations at the microscopic level are due to 
the thermal movement of adatoms surrounding the island 
constantly attaching to its edge and coming off from it. 

To diagonalize the free energy one rewrites the inte- 
grand in Fourier form 

F[{r n };t] = 27rA^ (An-n + mnB m - n ) r„(i)r^(t) , 

(7) 

where = J r(#) exp[ike]de and similarly for A), and 
-Bfe . The Fourier modes arc coupled due to the anisotropy, 
which is contained in A and B as we shall see shortly. 
Here n = is the expansion-contraction mode; n = 1, 
which we called the dipole mode in the Introduction, is 
related to the Brownian, diffusive motion of the island; 
n = 2 is a quadrupolar distortion, i.e. an elongated shape 
with two maxima and two minima in perpendicular di- 
rections; and so on. The Fourier components have her- 
mitian properties since A{9) and B(e), the factors associ- 
ated with the equilibrium island shape, are real functions; 
hence, = A*. _B_, = £?„*, and r_,- = rf. 
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FIG. 3: (Color online) Extension of Fig. 2, showing instanta- 
neous island shape, for analysis of fluctuations. For a particu- 
lar azimuthal direction, 8, the deviation from the equilibrium 
island shape, r, its derivative with respect to 9, r, the unit 
vector normal to the instantaneous shape, n, and the corre- 
sponding angle, ip, are all time dependent. 

The free energy of Eq. Q can readily be cast into 
matrix form: 

F[r; t] = 27rAr t (A + MBN) r , (8) 

where r is a vector containing the Fourier components of 
the instantaneous island shape, A and B are hermitian 
matrices, [A] m< „ = A TO _„, [B] m< „ = B m ^ n , and M = N 
are diagonal matrices with the wavenumbers along the 
diagonal. 

As in practice there are only a finite number of atoms 
on the edge of the island, we discretize the problem. If 
the number of atoms on the edge is 2N, there are as 
many modes in the system; as we will see in Sec. IV, to 
analyze the right number of modes is crucial to the prob- 
lem. Now, if r contains the Fourier components from 
— N + 1 through N, the Fourier transform is discrete and 
rfc = Y^,j=-N+i r j ex P[* kjn/N], where r| is the devia- 
tion from the equilibrium shape in the 9 — jn/N direc- 
tion. 

The Ak and the Bk can be obtained similarly; and A 
and B arc finite cyclic hermitian matrices, meaning that 
their diagonal elements are the same. They also reflect 
the symmetry of the equilibrium shape as e.g. in our case 
due to the six-fold symmetry the principal diagonal is 
filled with Aq, the 6th to the right with A-§, etc. As 
M and N are the same diagonal matrices, the MBN 
product keeps the hermitian property. 

In the isotropic case (when the equilibrium shape is 
circular), A(6) = and B(0) = 1/2 for all 9. After 
the Fourier transformation this gives A = (zero ma- 
trix) and B = (1/2)1 (diagonal matrix). The anisotropy 
comes into play when the equilibrium shape is not circu- 
lar, so that A(9) and B(9) are not constants and their 
higher order Fourier components fill the (off-) diagonals. 
These off-diagonals couple the Fourier modes. 



Due to hcrmiticity the above matrix form is diagonal- 
izable 

F[{h n };t] = 2TT\J2^nhnK, (9) 

n 

and the eigenvalues A„ of the A + MBN matrix are all 
real. As we see shortly (in Eq. 1)111 \) these eigenvalues 
are related to the strengths of the h n eigenmodes, which 
at every time instant are just the transforms of the r n (t) 
Fourier modes of the instantaneous island shape. Again 
due to hermiticity, there is a unitary matrix U which 
transforms Eq. JSJ into Eq. 10 and gives the linear rela- 
tionship between r and h: r = Uh, where the vector h 
contains the h n as its elements. 

This decomposition of the free energy into eigenmodes 
in Eq. © facilitates the calculation of the Lagrange mul- 
tiplier A. In equilibrium, according to the equipartition 
theorem, the ensemble average of each mode, represent- 
ing a degree of freedom, must have the same Boltzmann 
energy: 

E„ 

2ttAA„<|/i„| 2 ) = X -k B T. (10) 

A„ and (|/i n | 2 ) can be determined from the equilibrium 
island shape and the fluctuating island perimeter, respec- 
tively. Here E n must be a constant in n, the modes, as 
the temperature and the chemical potential, A, are fixed 
macroscopic parameters of the island. From this equa- 
tion one can determine the same A, in principle, from 
any mode. Thus, either experimentally observing island 
fluctuations or using Monte Carlo simulations one can 
determine E n . which in turn provides A. This A was the 
missing parameter to determine absolute line tensions, 
and plugging it back into Eq. J2J) , we get the line tension 
in all azimuthal directions. 



III. KINETIC MONTE CARLO 

The scarcity of extensive experimental data leads us 
to use Monte Carlo methods to simulate the behavior 
of the system. However, use of numerical rather than 
experimental data for testing of formal ideas has many 
advantages in any case. Most obviously, in numerical 
experiments one can obtain far greater control, with no 
worries about anomalous behavior due to unsuspected 
stray contaminants. Typically one can generate much 
more data. In the present experiment, we do not need 
to worry about the scan rate of the probe; our lattice 
configurations are instantaneous snapshots. Another ad- 
vantage of computer simulations will also become clear 
in the next section: it allows us to analyze correlation 
times. 

Since our original motivation was to simulate the re- 
laxation of a Pb crystallite with a (111) facet, we place 
a nanoscale island on a triangular lattice. We surround 
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TABLE I: Tabulation of some of the energy barriers used 
in KMC simulations of Pb(lll). The energies in columns 2 
and 3 were computed by M. Haftel using SEAM21 with glue 
potentials^ In the last column are the energies used in the 
simulations. (The edge-diffusion energy barrier in the Break- 
three scheme is closer to the corner rounding barrier of SEAM 
than to the actual straight-edge diffusion barrier of SEAM 
(see text); the latter has a much lower barrier: 108 meV; 
unfortunately, the bond-counting method does not distinguish 
between these two.) 



it by a non-permeable circular container of radius R c to 
let the system reach its thermodynamic equilibrium, in 
order to measure its equilibrium fluctuations Mdl Thus, 
this geometry corresponds to an island placed on top of a 
facet of a crystallite (with an infinite Ehrlich-Schwoebel 
barrier). Note that by adjusting the permeability one 
can tune the overall decay rate of the island, which in 
this paper we fix at zero. 

Since the temperature of the systems of interest is low 
compared to the energy barriers of adatomic hopping, 
we have chosen to use the Bortz-Kalos-Lebovitz (BKL) 
continuous-time MC algorithm^ as it is best suited to 
low temperature systems and as its rejection- free method 
allows us to greatly improve the efficiency of the simu- 
lations compared to traditional Metropolis algorithms. 
The typical temperatures are around T c /4 or less for the 
two-dimensional lattice gas of adatoms on the surface. 
Using the n-fold way method to keep track of the avail- 
able MC moves, we could improve the efficiency even 
further. Because of the small number of energy barriers, 
the n-fold way approach (5-fold) is superior to the binary 
tree implementation of the BKL algorithm^ in this case. 

Since we are not interested in all the details of this sur- 
face in the simulations, but only try to capture the main 
mechanisms, we do not take into account the ABC stack- 
ing structure of the fee lattice of Pb(lll). Hence, the top 
layer constitutes a triangular lattice with perfect six-fold 
symmetry. Furthermore, we assume that adatoms can 
only hop to nearest neighbor sites, and that the energy 
barrier the adatom has to overcome is determined by the 
occupation of the 8 sites surrounding, as nearest neigh- 
bors, the 2 sites involved in the hopping process. 

The energy barriers for hopping rates are mainly based 
on the semiempirical embedded atom method (SEAM)21 
using Ercolessi's glue potentials^ for the Pb(lll) surface 
to derive the bond-counting energy barriers^ actually 
used in the simulations (see Table HJ) . We use variants of 
simple bond-counting. In what we term the break-three 



scheme, we count only bonds that have been broken with 
the three sites "behind" the move. If we also include 
bond-breaking with the 2 sites to the left and the right 
of the move (denoted side sites), we have a break-five 
scheme. The 3 sites in front of the move do not affect 
the energy barrier in our simulations. The break-three 
and break-five schemes both satisfy detailed balance in 
a straightforward fashion. They both should give the 
same results for static parameters, since they are both 
nearest-neighbor schemes. Comparison runs using this 
feature provided one test (among many) of our program. 
However, the kinetics obviously differ because the energy 
barriers tend to be higher for the break-five scheme, slow- 
ing the kinetics significantly. Since the energy barriers in 
the break-three case are closer to the calculated barriers, 
we choose to use this scheme in our simulations. 

In Table [I] we list a few energy barriers calculated 
by the above-mentioned SEAM and the corresponding 
break-three-scheme barriers. Surface diffusion is when 
an adatom has no other adatoms in the surrounding 8 
sites in SEAM calculations; in the break-three scheme 
this barrier also applies to all cases in which any of the 
side sites or the sites in front are occupied. This is the 
reason why the attachment barrier is the same as the 
surface diffusion barrier. Edge diffusion is the case in 
which a side site is occupied, as is its nearest neighbor 
"behind" . The energy barrier associated with this pro- 
cess is 237 meV whereas if there is a nearest in the front 
as well so that the adatom rolls along three others on 
one side the third in the front seems to assist the hop a 
great deal (at least according to the SEAM data) as the 
barrier is 108meV. The break-three scheme has the same 
barrier for these two processes and also for any other in 
which only one bond is broken. The break-two-bonds 
barrier corresponds to a hop with 2 occupied sites in the 
back, Break-3-bonds is when all 3 sites are occupied in 
the back. The very high "Out" energy barrier assures 
that adatoms cannot escape the container; i.e., the per- 
meability is zero. 

The basic parameters of the surface investigated and 
the KMC simulations are the following: The nearest- 
neighbor spacing a\ on the Pb(lll) surface is 3.50A. 
The typical island diameter 2R is 40ai to 80ai, while the 
container diameter R c ranges from 12.5% to 300% larger 
than the island. We examined temperatures 250K, 300K, 
350K, and 400K. In each MC snapshot of the island, we 
measure the island radius from the instantaneous center 
of mass in N "equiangular" directions where N = 360 if 
not indicated otherwise. 

We start the simulations from a nearly circular shaped 
configuration and let it relax to equilibrium, starting 
the MC measurement of the fluctuation and shape after 
the longest wavelength mode has passed its correlation 
time. Especially at the lower temperatures, the typical 
equilibration times are very long, consistent with other 
reportsiSLSi Ensemble averages are taken from 100 to 
3000 different runs starting from the same initial con- 
figuration, but with different random-number seeds. In 
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FIG. 4: (Color online) Eigenmodes (open squares) and 
Fourier modes (solid circles) at T=250K, R=20a±, R c =4Qa\. 
E n is measured in atomic spacing units. 



each run, after equilibration, we get statistically inde- 
pendent fluctuations at time intervals again determined 
by the relaxation time of the longest wavelength mode. 
We take such independent "snapshots" of the islands 5 
to 200 times in each run, so that we typically have 10,000 
to 70,000 islands over which to average. 



IV. CHEMICAL POTENTIAL, LINE TENSION, 
LINE STIFFNESS 



As described in detail in Sec. II, the energetic parame- 
ters of the island edge are determined by the island shape 
and its edge fluctuation. The Wulff construction pro- 
vides the relationship between the relative line tension in 
the azimuthal directions on the crystal surface and the 
equilibrium island shape, and the information from the 
fluctuations {\h n ^) of each mode in Eq. jTDJ> gives the 
chemical potential A that makes the Wulff construction 
absolute in Eq. 

From our KMC simulations we determine E n of Eq. 
ffTUft: it is depicted in Fig. E|for T=250K and R=2Qa x 
island radius. Since the perimeter is about 120ai, we use 
120 points to describe the circumference out of the 360 
available. 

We calculate E n both using the transformation to the 
eigenmodes taking into account the anisotropy, and also 
pretending the islands were isotropic. In this latter case 
the h n = r„ are simply the Fourier modes, and A„ = 
(l/2)n 2 . 

The Fourier modes and the eigenmodes are nearly in- 
distinguishable (cf. Fig. except for long wavelengths 
(small wavenumber n) . The "fluctuations" in the Fourier 
modes for longer wavelengths are the only signature of 
the anisotropy. This effect is smoothed out by transform- 
ing to the eigenmodes. (Only the longest [n=2] wave- 
length mode seems to stand out after the transforma- 
tion. We arc still puzzled that the longest wavelength is 
so special and does not couple to the shorter wavelength 



modes. We give a more detailed account of the analy- 
sis of the coupling elsewhere The chemical potential 
can be determined from either the Fourier modes or the 
eigenmodes using their plateau regions (between n = 4 
and n = 40); however, since only the longest- wavelength 
modes (if any) are measurable due the poor temporal res- 
olution of present-day experimental apparatus (see the 
results for correlation times in the next section), the 
transformed modes serve better for determining chem- 
ical potentials. Finding E n thus from these intermediate 
wavenumbers, gives E n =0.275af and through Eq. I|10|) 
and Eq. J2J /? =34.1 meV/A for the line tension for the 
high-symmetry direction. This value reasonably approx- 
imates the experimentally obtained ones for Pb(lll) at 
T=393 K: (3 1A = 27.9 meV/A and /3 1B = 26.5 mcV/A for 
A- and B-type steps, respectively, considering the crude 
approximation of bond-counting mentioned earlier^i In 
our simulations the two directions corresponding to the 
two different types of steps are intrinsically equivalent be- 
cause we assume six- fold symmetry as the available values 
for energy barriers that we use in our KMC do not dis- 
tinguish between the A- and B-dircctions, as mentioned 
in Sec. III. 

At higher temperatures the Fourier modes deviate less 
from the eigenmodes for longer wavelengths, as expected 
since the equilibrium shape is more nearly circular and 
less affected by the underlying anisotropy. 

In earlier workA&S^ in which experimental data were 
used as an input of similar calculations, there is a sum 
over the modes, but because those modes are buried in a 
sum in the generalized cquipartition theorem, one cannot 
see whether they are the modes which satisfy, at least to 
a certain extent, the equipartition theorem. In those ex- 
perimental data the correlation times of the modes are 
not known and the effect of the finite temporal resolu- 
tion of the experiment may also interfere with the fluc- 
tuations which should in principle be determined from 
"snapshots" , i.e. fast scanned images — fast at least com- 
pared to the correlation times of the modes used in such 
calculations. We shall elaborate on this in the next sec- 
tion. 

For the same temperature we do the same measure- 
ment as above, but monitored N = 180 points on the 
perimeter instead of N = 360. The Fourier modes are de- 
picted in Fig. [3 There are approximately 120 atoms on 
the perimeter, but since we cannot divide the 180 perime- 
ter points into 120 equiangular ones to do a Fourier trans- 
forms, we use A^=90 or ^=180 perimeter points as an 
approximation and observe how the plateau changes from 
what we saw in Fig. 0] The comparison of these two 
plots from MC simulations might help analyzing exper- 
imental data with limited spatial resolution as well, as 
it shows how the Fourier modes behave in case of un- 
dersampling (A^=90) and oversampling (^=180). The 
undersampled modes give higher values for E n than ex- 
pected for modes \n\ > 25, as if those modes took over 
the energy of the modes that are missing from spectrum 
(namely 45 < \n\ < 90). In oversampling there is not 
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FIG. 5: Fourier modes for N=90 points (open squares) and 
iV=180 points (solid circles) on the perimeter at T=250K, 
and R—20ai, R c =40ai. E n is measured in atomic spacing 
units. 



enough energy for all modes in the sampling, so they go 
below the expected value of E n . This is the simple rea- 
son for the peculiar shape of the two curves in Fig.[S]The 
value of E„ can still be determined quite accurately by 
using the value at which the two curves start to separate 
at \n\ ~ 10, providing basically the same result for E n as 
before. 

From the equilibrium island shape using the Wulff con- 
struction, we have determined the relative line tensions 
and stiffnesses in the azimuthal directions on the (111) 
surface (see Fig. |SJ). The equilibrium shape is more and 
more "faceted" in the six main directions as we can ex- 
pect, the stiffness is about 3 times bigger in the direc- 
tion of the "facet" than in the direction of the corner for 
T=350 K, whereas this factor is about 20 for T=250 K, so 
that it spikes out much more, but still has only a smaller 
effect on the spectrum as shown before. 



AUTOCORRELATIONS, CORRELATION 
TIMES, KINETICS 



Inspecting the autocorrelation functions in Fourier 
space, we find that the longest wavelengths have sur- 
prisingly long correlation times (in CPU time). Thus, 
for our fairly large (at least for the computer) systems, 
it is hard to reach full equilibration needed to make the 
desired MC measurements. Most surprising is that the 
relaxation time of the longest-wavelength modes is 10 to 
100 times longer than the relaxation time of the islands 
to their equilibrium shape. Hence, estimating the ther- 
malization time from just the shape relaxation may be 
very misleading, possibly giving problematic results not 
characteristic of equilibrium. Such behavior may include 
illusory strong-mode coupling, or stronger fluctuations in 
autocorrelation functions even in case of good statistics. 



The temporal correlations can be characterized by 

G{t) = ([r{tv)-r{t a +t)] 2 )<xt 2 P. (11) 

Since we measure correlations in equilibrium, to must 
be greater than the thermalization time of the system. 
Here r(t) is the fluctuation from the equilibrium shape, 
as before, and depends on the angle, 6, and time. The 
growth exponenti%3 characterizes the temporal behavior 
of the fluctuations. The average is taken over angles and 
an ensemble as well. 

The typical behavior of the correlation function is that 
the exponent /? remains at 1 for very short times^ typ- 
ical of the ballistic behavior of diffusion at very short 
times, and then crosses over to a value which character- 
izes the rate-limiting kinetics driving the fluctuations of 
the island edge; eventually it crosses over to zero as the 
correlation function saturates due to the finite size of the 
system. 

Pure rate-limiting kinetics have been thoroughly 
investigated>ii*i^2iS4^4 In these well-defined cases, /3 
can take the values 1/4 for attachment-detachment ki- 
netics, 1/6 for surface diffusion, and 1/8 for step-edge 
diffusion, where the last mechanism gives a very "slow" 
dynamics. There can be crossover regimes between these 
pure cases, leading to values of (3 between the quoted 
values, and certain geometries can also effect the value 
of (3. One should also see crossovers as length scales 

To investigate the length-scale dependence of the cor- 
relation function, it is more appropriate to use the cor- 




FIG. 6: (Color online) Polar plot of the equilibrium island 
shape R(0) (outer dots), the relative line tension /3(ip) (inner 
dots) and the relative line stiffness f3(ip) (the innermost curve) 
in arbitrary units at T =350K, R=2Qa\, R c —40ai. (note the 
difference between 9 and ip) 
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FIG. 7: (Color online) Correlation function G n (t) of the 
Fourier modes for n=2,3,...,10 from top to bottom. T =400K, 
H=20oi, R c =40ai. 

relation function in Fourier space: 

G n (t) = (|r„(t ) - r„(t + t)\ 2 ) (12) 
= C n (1 - exp (- \t\ /r n )) , (13) 

where the C n are twice the amplitudes of the fluctuations 
of the modes analyzed in the previous section, and the r„ 
are their correlation times. The wavenumber dependence 
of r n is known to have an intimate relationship with the 
exponent in G(t), namely 

r„~n- z , (14) 

where z is the dynamic exponent, and the scaling rela- 
tionship between z and /3 here is z = (l/2)/[3w& The cor- 
relation time increases with increasing wavelength, with 
the scaling exponent z. For larger exponent z, the corre- 
lation times grow more rapidly, so that for longer wave- 
lengths the correlations and the dynamics in general can 
slow down very "quickly" . 

Here we pay particular attention to the longest wave- 
length and its corresponding correlation time T2, which 
makes the largest contribution to the fluctuations and re- 
laxes the most slowly. From the wavelength dependence 
of the correlation time, we also determine the dynamic 
exponent and the rate limiting kinetics. 

In the KMC simulations for T=400K and R=20a 1 , the 
longest-wavelength mode, n=2, is 120ai or 420A long. 
From Fig. its relaxation time is T2 = 5.5 x 10 7 MCS 
(Monte Carlo steps). To give a crude estimate for T2 in 
real time, we consider the hopping rate 

v = v D exp[-/3E b ] (15) 

to be the product of the attempt frequency, which we 
identify with the Debye frequency of Pb: vq = 1.83xl0 12 
Hz^ and the Boltzmann factor of the energy barrier of 
a particular hop. Hence, a MCS in this Monte Carlo 
simulation is equivalent to a 1/vd time increment in real 



time; thus, the relaxation time in this particular case is 
t-2 = 0.030 msec. 

As expected, these correlation times change dramati- 
cally with temperature as the underlying physical phe- 
nomena are activated. For R=20a\ at T=350 K, T2 = 
2.0 x 10 8 MCS or 0.11 msec, which means 4 times 
longer relaxation compared to 400 K, while for 300 K 
T2 = 7.1 x 10 s MCS or 0.39 msec, which represents slow- 
ing by another factor of 4. 

We have to mention here that we did two sets of 
kMC simulations. The simulations described in this pa- 
per satisfy detailed balance and sample the canonical 
distribution, while in the other set we choose the en- 
ergy barriers to be those calculated by SEAM, which 
cxplicitely violate detailed balance and energy conser- 
vation. Interestingly, the latter simulations tend to give 
closer agreement with experimentally measured correla- 
tion time data. We do a comparison of the results of 
these two sets of simulations^ and upcoming experimen- 
tal data^i elsewhere. 

The scaling of the relaxation time with (wave)length 
can be seen in Fig. [H] In the plotted wavenumber 
range, overall, r« basically behaves like z =4, sug- 
gesting that the mechanism driving the fluctuations is 
step-edge diffusion which is in agreement with previous 
observationsi^SiS 7 . 

Comparison of these length scales and their corre- 
sponding relaxation times with existing experimental ob- 
servations might give interesting physical insight. For 
example, in the experiment by Thurmer et al.^ a small 
Pb crystallite, whose top facet has a perimeter slightly 
larger than 1200 nm, relaxes at 383K to its equilibrium 
(or at least metastable) shape in 1-2 days after being 
quenched from a higher temperature. Since the crystal- 
lite is 30 times larger than the longest wavelengths in 
our KMC simulations, the relaxation time T re i ax of the 
longest wavelength mode of the perimeter of the top- 
most island on the crystallite is 8.1 x 10 5 times longer if 
the rate-limiting kinetics is step-edge diffusion (though, 
of course, attachment-detachment and terrace diffusion 
could be present but NOT rate- limiting). Specifically, 
the value of r re i ax is 24.3 sec based on our KMC data at 
T=400K. From these data it seems that the island fluctu- 
ation is much faster than the decay of the 3D structure; 
thus there is no direct relationship between the fluctua- 
tion and the decay. 

The above arguments lead to a general view of the 
evolution of surface structures. For Pb in the temper- 
ature range 350 K-400 K, one observes the slow devel- 
opment and relaxation of fluctuations at the /im scale 
in experiments. Assuming that the rate- limiting kinetics 
retain the same z =4 range for even longer wavelengths, 
structures of 10 (im size — step edges, islands, etc. - 
take days to years to change due to the large dynamic 
exponent, z, so in effect they look frozen under labora- 
tory conditions. This is the reason why these monolayer 
structures do not show any large-scale changes while on 
a shorter scale they can be very active. The structural 
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FIG. 8: (Color online) Correlation time r n vs. wave num- 
ber n on a log-log scale. The MC data (black circles) show 
z = 4 which corresponds to step-edge diffusion. The dashed 
lines represent z=4, 3, and 2 (from top to bottom) dynamic 
exponents. T =400K, 7?=20ai, R c =AQa%. 



changes in the 3D crystallites are even slower, they are 
even more stable. 

Lowering the temperature makes the length scales — 
at which evolution or relaxation can be observed — ex- 
ponentially shorter, which is readily understandable if 
one looks at the converse of the above arguments. The 
length scales with time like I ~ r" 1 / 4 (for z = 4) whereas 
r scales like \/v in Eq. (|15|l . Thus, the length scales 
with temperature like I ~ exp {E^/AksT]. This basically 
means that given the temperature and the time scale of 
observation, one can calculate an effective length scale, 
I/ e ff on which the structures on a surface are in equilib- 
rium with their surroundings and actively changing on 
the time scale of e.g. STM measurements. Having L e s 
and the energy barrier Ef, at a certain temperature, one 
can also make at least rough estimates of the effective 
lengths at other temperatures using this scaling argu- 
ment. This picture is certainly a result of simplification, 
since there is a whole set of energy barriers in such a 
complex physical system as a crystal surface, and the 
various atomistic mechanisms governed by different bar- 
riers freeze out or get activated at different temperatures, 
depending on their corresponding energy barriers^ 

Comparison or extrapolation to other materials is pos- 
sible if the energy-barrier set is similar to that of Pb ( 1 1 1 ) . 
Then the Debye frequency sets the time scale while the 
energy barriers set the temperature scale, as one can 
readily deduce it from Eq. I|15|) . On the other hand, if the 
energy-barrier set is completely different, as for example 
for Si^S it gives rise to a different rate-limiting mecha- 
nism namely attachment- detachment for a wide range of 
temperatures, and such extrapolation is not possible, but 
a whole set of simulations should be done for the group 
of materials with this sort of barriers. 



VI. CONCLUSIONS 

In this paper we deduce energetic and kinetic param- 
eters of a particular metal surface below its roughening 
temperature. We use kinetic Monte Carlo simulations 
to mimic the fluctuations of large nanoscale islands on 
these smooth surfaces in order to determine equilibrium 
island shapes, anisotropic line tensions in the azimuthal 
directions of the surface, and the correlation times of the 
Fourier modes of the fluctuations. 

We derive an analytic expression for the chemical po- 
tential of the island edge from the equilibrium island 
shape and the associated capillary wave fluctuations 
around it. This chemical potential sets the scale for 
the anisotropic line tension, [the azimuthal dependence 
of] which is usually known only up to a multiplicative 
constant. To account for the anisotropy of the line ten- 
sion, this procedure contains a transformation from the 
Fourier modes of the island edge fluctuations to the true 
eigenmodes. However, detailed analysis of the Fourier 
and eigen modes of the fluctuations reveals that the dif- 
ference in their spectrum (Fig. [3} is unexpectedly small. 

The obtained line tensions — one of the most impor- 
tant physical parameters of steps on surfaces — are in the 
correct range compared to known experimental results, 
even in this simplistic model, with its rather small set of 
hopping-cnergy barriers in the KMC simulation. 

We have analyzed the effect of spatial sampling, which 
shows that the long wavelength modes are hardly af- 
fected by the undcrsampling (oversampling) — too low 
(too high) resolution of imaging, which means that there 
are fewer (more) sample points on the step edge in the 
image than actual atoms in the experiment — whereas 
moderately short and short modes change significantly. 

The analyses of the correlation times of the Fourier 
modes show that nanoscale objects fluctuate on the msec 
to /xsec time range at moderately high temperatures 
(400K) on Pb surfaces. Since the atomic processes are 
activated, this time scale changes dramatically with tem- 
perature. 

In closing, we comment on the equilibration time of 
step structures in Monte Carlo simulations. The full equi- 
libration of these structures is signalled by the correlation 
time of the longest wavelength mode, which can be very 
large (in CPU time) for the system sizes and temper- 
atures studied in this paper. To do correct MC mea- 
surements in equilibrium, one has to pass this time; oth- 
erwise, results for "equilibrium quantities" can be very 
misleading, as is well known from non-equilibrium sta- 
tistical mechanics. If one docs not look at correlation 
times of Fourier modes, very careful analysis is required 
to avoid such equilibration problems^ Recently, several 
papers have appeared concerning correlation functions, 
persistence, etc., of steps much longer than ours^ and 
sometimes even several of them (in studies of the inter- 
action between them) . They might well suffer from these 
problems since this equilibration time scales with system 
size as the fourth power, meaning that a system that is 
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twice as big needs an order of magnitude longer CPU 
time to be equilibrated. 

We also point out that the measurement of these fluc- 
tuations experimentally might be difficult because of the 
above-mentioned time scale of the fluctuations. One ei- 
ther has to use techniques with which snapshots of the 
surface can be taken^or, in direct visualization methods 
(like STM measurements) , the scan rate of the equipment 
must be faster than the fluctuations of a given wavelength 
of interest. Otherwise one measures the two ends of a 
wavelength at such a time separation that they are un- 
correlated, leading to difficulties of interpretation. The 
"speed" of the fluctuations can be tuned by changing the 
temperature, but one also has to take into account that 
lowering the temperature decreases the size of the fluc- 
tuations, rendering the measurement harder. 

Finally, the extrapolation of our results for nano- 
objects to mesoscale features makes possible comparisons 
of correlation times of modes of certain wavelengths, as 
well as of decay or of relaxation of larger structures to 
their equilibrium forms. This comparison reveals that 



the relaxation of 3D structures of the same lengths are 
slower than that of the simple step or monatomic high 
islands, due to additional mechanisms and phenomena 
like a possible Ehrlich-Schwoebel barrier, very low de- 
tachment rate and perhaps elasticity that affect the 3D 
structure. 
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